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MATRIX  ANALYSIS  OF  LONGITUDINAL  AND  TORSIONAL  VIBRATIONS 
IN  NONUNIFORM  MULTIBRANCH  BEAMS 


By  Robert  T.  Wingate 
Langley  Research  Center 

SUMMARY 

Since  longitudinal  modes  and  frequencies  provide  basic  data  for  dynamic  analyses 
of  arbitrary  beam-like  structures  and  since  closed-form  solutions  for  the  modes  are 
generally  not  feasible  to  obtain,  an  approximate  method  is  developed  for  computing 
the  natural  frequencies  and  the  corresponding  mode  shapes  for  a  variable- section, 
unconstrained  multibranch  beam. 

A  lumped  mass  analogy  employing  influence  coefficients  is  used  to  represent  the 
beam.  The  simultaneous  equations  of  motion  for  the  lumped  mass  system  are  derived  in 
matrix  form  and  algebraically  manipulated  to  yield  a  classical  eigenvalue  equation  solv¬ 
able  by  standard  procedures.  The  orthogonality  relationship  of  the  natural  modes  is 
derived  and  used  to  form  the  basis  of  an  orthogonal  sweeping  process  for  determination 
of  modes  above  the  fundamental. 

Numerical  examples  including  an  application  to  a  solid-fuel  launch  system  are  pre¬ 
sented.  Also,  a  detailed  discussion  is  devoted  to  the  theoretical  verifications  of  the 
approximate  modes  and  frequencies. 


INTRODUCTION 

The  capability  to  calculate  the  longitudinal  natural  frequencies  and  associated  mode 
shapes  of  beam-like  structures,  such  as  solid-fuel  launch-vehicle  systems,  railroad 
trains,  machine  shafts,  and  piles,  is  of  basic  importance.  The  natural  frequencies  pro¬ 
vide  necessary  design  information  and  the  modal  functions,  because  of  their  special 
orthogonal  properties,  are  ideal  for  use  in  "series"  solutions  to  the  differential  equations 
of  motion.  Such  "series"  or  "normal  mode"  solutions  provide  a  powerful  mathematical 
tool  for  the  analysis  of  transient  or  steady-state  response  of  a  structure  to  disturbing 
forces. 

For  purposes  of  calculating  the  frequencies  and  mode  shapes,  slender  structures  of 
the  type  mentioned  are  usually  represented  as  a  simple  beam  which  may  have  a  varying, 
even  discontinuous,  cross  section.  However,  for  structures  carrying  elastic  appendages 


(e.g.,  a  launch  vehicle  with  a  flexible  payload  within  a  flexible  heat  shield)  it  is  sometimes 
necessary  to  represent  the  structure  as  a  branched  beam  -  that  is,  a  main  beam  with 
subsidiary  beams  attached  at  points  along  its  axis. 

Exact  solutions  for  the  longitudinal  frequencies  of  even  a  simple  beam  are  not  fea¬ 
sible  except  for  very  special  cases  of  the  cross-sectional  variation  (the  solution  for  a 
constant  cross  section  is,  of  course,  well  known).  In  reference  1  Bessel  function  solu¬ 
tions  were  obtained  for  a  limited  number  of  variable-cross-section  beams  with  free-free 
ends  and  fixed-free  ends.  Also,  in  reference  2  Bessel  function  solutions  were  obtained 
for  a  fixed-free  beam  with  a  linearly  varying  cross  section.  In  practical  work,  however, 
recourse  is  generally  to  approximate  solutions. 

The  purpose  of  this  paper  is  to  present  an  approximate  method  for  computing  the 
longitudinal  natural  frequencies  and  mode  shapes  for  a  general  variable- cross- section, 
multibranch  beam.  The  method  is  derived  for  free  or  unconstrained  beams;  however,  the 
extension  of  this  procedure  to  other  boundary  conditions  is  straightforward.  Through  a 
change  of  variables  the  derived  equations  are  also  applicable  to  torsional  problems  in 
which  the  branch  beams  are  everywhere  concentric  with  the  main  beam. 

Available  approximate  methods  for  the  calculation  of  frequencies  and  modes  of 
structures  are  voluminous,  but  literature  on  their  application  to  longitxidinal  vibrations 
of  multibranch  structures  is  scant.  The  longitudinal  vibrations  of  beams  have  been 
treated  in  several  papers  (e.g.,  ref.  3)  by  the  so-called  "stiffness”  matrix  approach  in 
which  the  beam  is  idealized  as  a  number  of  point  masses  connected  by  springs.  The 
simultaneous  equations  of  motion  for  the  point  masses,  when  written  in  matrix  form,  yield 
a  classical  eigenvalue  problem  which  is  solvable  by  standard  methods.  This  approach  is 
efficient  when  it  works  but  suffers  from  the  disadvantage  that  the  lower  frequencies,  which 
are  generally  of  prime  interest,  are  associated  with  the  subdominant  eigenvalues.  The 
author’s  experience,  as  well  as  that  of  others  (ref.  4),  has  shown  that  matrix  eigenvalue 
problems  break  down  in  the  computation  of  the  subdominant  eigenvalues  whenever  they 
become  small  in  comparison  with  the  dominant  eigenvalue. 

The  spring-mass  idealization  of  a  structure  has  also  been  utilized  (e.g.,  ref.  5)  to 
analyze  longitudinal  vibrations  by  a, method  associated  with  Holzer  (ref.  6).  This  approach 
has  an  advantage  for  structures  with  many  discontinuities  in  the  stiffness  and  mass  in 
that,  through  a  recurrence  equation,  a  very  large  number  of  mass  points  may  be  used 
without  storing  a  large  matrix  in  the  computer.  A  fine  definition  of  the  mode  shape  can  . 
thus  be  obtained.  It  has  the  disadvantage  that  it  is  prone  to  roundoff  error  and  other  loss 
of  numerical  significance  especially  in  the  higher  modes.  Also,  close  attention  is 
required  to  avoid  missing  modes  in  the  trial  and  error  search  for  frequencies,  and  a  large 
number  of  iterations  are  required  to  cover  the  frequency  spectrum  of  all  the  modes. 
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The  method  presented  in  this  paper  is  a  "flexibility"  or  influence  coefficient  matrix 
approach  for  the  analysis  of  natural  longitudinal  or  torsional  vibrations  of  a  multibranch 
beam.  This  method,  which  was  originally  developed  in  reference  7,  is  an  adaptation  of  a 
procedure  described  in  reference  8  for  computing  the  lateral  bending  modes  and  frequen¬ 
cies  of  a  single  beam.  It  is  attractive  because  of  simplicity  in  derivation  and  its  adapta¬ 
bility  to  digital  programing.  The  major  advantage,  however,  is  that  the  lower  frequencies 
are  associated  with  the  dominant  eigenvalues  and  potential  problems  of  accuracy  loss  in 
the  lower  modes  and  other  losses  in  numerical  significance  are  virtually  eliminated.  The 
development  is  completely  general  and  requires  only  basic  parameter  changes  from  one 
application  to  another. 

By  proceeding  from  basic  physical  assumptions,  simple  matrix  algebra  is  used  to 
obtain  the  characteristic  equation  in  matrix  form.  Iterative  procedures  for  obtaining  the 
frequencies  and  modes  are  then  discussed.  Numerical  examples  are  given  for  an  ideal¬ 
ized  beam  and  a  typical  solid-fuel  launch  vehicle  to  illustrate  the  formulation  and  solution 
of  the  derived  matrix  equations.  Also,  a  detailed  discussion  is  devoted  to  the  accuracy  of 
the  method  in  comparison  with  exact  solutions  and  solutions  by  the  "stiffness"  matrix 
approach. 


SYMBOLS 


A  stressed  cross-sectional  area,  inches^  (meters^) 

A(j)  stressed  cross-sectional  area  of  jth  branch,  a  function  of  x,  inches^ 

(meters2) 

E  Young's  modulus,  pounds  force/inch^  (newtons/meters^) 

E(j)  Young's  modulus  of  jth  branch,  a  function  of  x,  pounds  force/inch2 

(newtons/ meter2) 


^r(j)’^s(k) 


axial  force  at 


r(j) 


Xg(k),  pounds  force 


(newtons) 


Gr(j)  influence  coefficient,  deflection  relative  to  of  rth  mass,  of  jth 

branch,  due  to  a  unit  force  at  sth  mass  of  kth  branch,  inches/pound  force 
(meter  s/new  ton) 


j  index  denoting  jth  branch 

k  index  denoting  kth  branch 
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L 

m 


^s(k) 


M 


overall  length  of  main  or  0th  branch,  inches  (meters) 


mass  per  unit  length,  pounds- seconds^  /newlor^-_se_coi^^^^^ 

inch^  \  meters^ 


sth  mass  of  kth  branch  located  at  x 
newtons- seconds^ 


s(k)' 


pounds- seconds^ 


inch 


total  mass, 


meter 

pounds- seconds^  /newtons- seconds^' 


inch 


meter 


n  index  denoting  total  number  of  lumped  masses  in  a  system 


P 

q 

t 

u 


v(k) 

X 

''r(i)’^s(k) 

V 


index,  corresponds  to  pth  natural  mode 
index,  corresponds  to  qth  natural  mode 
time,  sec 

displacement  from  unstrained  position,  function  of  x  and  t,  inches 
(meters) 

index  denoting  total  number  lumped  masses  on  kth  branch 
longitudinal  coordinate,  inches  (meters) 

distance  of  r(j)th,  s(k)th  lumped  mass  from  origin,  inches  (meters) 
exponent  for  parameter  study 


X  eigenvalue 

^  longitudinal  coordinate,  inches  (meters) 


‘^r(j)’‘^s(k) 


displacement  amplitude,  function  of  x,  denotes  deflection  from  unstrained 
position  of  r(j)th,  s(k)th  discrete  mass,  inches  (meters) 


CO  circular  frequency  of  simple  harmonic  vibration,  radians/ second 
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Subscripts: 


a  ath  row  of  a  matrix 

b  bth  column  of  a  matrix 

r(j)  rth  mass  of  jth  branch 

s(k)  sth  mass  of  kth  branch 

max  maximum 

ref  reference 

Matrix  notation: 


diagonal  matrix 


L  ^ 


square  or  rectangular  matrix 
T 

transpose  of  matrix 


1 

1 

1 


row  matrix 

column  matrix 

identity  matrix 

square  matrix  of  all  unit  elements 
row  matrix  of  all  unit  elements 

column  matrix  of  all  unit  elements 


ANALYSIS 


A  lumped  mass  analogy  of  a  beam  is  developed.  This  analogy  is  then  extended  to 
obtain  a  lumped  mass  representation  of  a  general  unconstrained  multibranch  beam. 
Through  an  influence  coefficient  approach,  the  longitudinal  equations  of  motion  for  the 
lumped  mass  system  are  derived  in  terms  of  one -dimensional  beam  theory  and  reduced 
to  a  matrix  eigenvalue  problem  which  determines  the  natural  frequencies  and  the  asso¬ 
ciated  modal  vectors.  A  procedure  for  computing  the  influence  coefficients  is  given,  the 
orthogonal  properties  of  the  modal  vectors  are  demonstrated,  and  an  iterative  method  is 
presented  for  solving  the  eigenvalue  problem.  Since  these  equations  are  based  on  one¬ 
dimensional  beam  theory,  they  should  be  applied  with  reservation  to  structures  which 
depart  from  this  concept  (e.g.,  thin-wall  cylindrical  shells  of  large  diameter,  shells  con¬ 
taining  a  sloshing  fluid,  etc.). 

Derivation  of  the  Lumped  Mass  Equations  for  a  Single  Beam 

The  so-called  ^’lumped  mass’’  analogy  in  which  a  structure  is  represented  as  a 
series  of  lumped  masses  connected  by  massless  springs  offers  an  effective  simple 
approach  to  vibration  analyses.  The  transition  from  a  continuous  system  to  a  lumped 
system,  however,  has  generally  been  somewhat  of  an  arbitrary  process  which  depended 
on  the  judgment  of  the  analyst.  A  formal  derivation  of  the  lumped  mass  approach  is 
obtained  by  the  direct  collocation  method  as  shown  herein.  This  process  involves  satis¬ 
fying  the  integral  equations  of  motion  at  a  finite  number  of  selected  points. 

Consider  a  general  nonuniform  beam  (sketch  1)  which  is  fixed  at  an  arbitrary  point 
Xj  with  a  concentrated  axial  unit  load  applied  at  X  =  The  resulting  deflection  at  any 
point  X  due  to  the  unit  load  at  |  is  designated  as  the  Green’s  function  or  influence 
function  G(x,^). 


Sketch  1 

By  utilizing  this  influence  function,  the  deflection  u(x,t)  at  any  point  x  due  to  a 
distributed  axial  force  F(^,t),  which  includes  inertial  forces,  is  given  by  superposition  as 


If  the  fixity  constraint  at  xj  is  released  and  the  beam  is  allowed  to  translate  as  a  free 
body,  equation  (1)  becomes 

u(x,t)  =  u^Xf,t)  +  G(x,^)  F(^,t)  d|  (2) 

where  u^Xf,t)  is  the  rigid  body  displacement  of  the  influence  function  reference  point  xf. 
Assume  that  the  imconstrained  or  free-free  beam  vibrates  in  simple  harmonic  motion  of 
frequency  w  and  amplitude  (p{x);  then, 

u(x,t)  =  (p{x)  e^^^*  (3) 

where  i  =  In  the  absence  of  external  forces  the  axial  force  is  caused  only  by  the 

inertia  loading  of  the  distributed  mass  m(x)  and  is  given  as 

F(x,t)  =  -m(x)  u(x,t)  =  co^m(x)  (p{x)  e^^^  (4) 

where  a  dot  denotes  differentiation  with  respect  to  time.  Substitutii^  equations  (3)  and 
(4)  into  equation  (2)  and  canceling  e^^^  gives 

(p(x)  =  <iP(Xf)  +  G(x,|)  m(^)  cp{^)  d^  (5) 

The  beam  length  from  x  =  0  to  x  =  L  is  divided  into  v  intervals  (sketch  2) 
which  are  arbitrarily  spaced  except  that  an  interval  marker  is  assumed  to  coincide 
with  any  discontinuity  in  axial  stiffness.  The  interval  boundaries  are  denoted  by  Xg 
where  s  *=  1  to  v  +  1. 


Xl  X2  Xg  Xy  X^^^ 

Sketch  2 

By  introducing  the  function  H(x,  =  G(x,^)  equation  (5)  can  be  written  as 


V 

cp(x)  =  cp^Xf'j  +  Y  ^  H(x,|)  m(|)  d| 


The  structure  is  continuous  in  the  neighborhood  of  Ig  so  that  the  function  H(x,  |)  can  be 
expanded  into  a  Taylor  series  about  ,^g,  where  Xg  <  |g  <  Xg^j^,  to  give 
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(7) 


H(x,^)  =  H(x,4g)+  h’(x,?s)(|  -  Is)  +  OA|s^ 

where  the  prime  denotes  d/d|  and  OA|g^  denotes  terms  of  order  -  |g)^  and 
higher.  Neglecting  higher  order  terms  and  substituting  equation  (7)  into  equation  (6) 
yields  the  approximate  relationship 


<p(x)  =  <p(xf)  +  y  H(x,|g)  m(|)  d|  +  H’(x,|g)  m(|)  (|  -  |g)d| 

s=l  ^s 


(8) 


In  equation  (8)  if  |  is  chosen  to  be  the  mass  center  of  gravity  between  Xg  and  Xg^j, 
s 

then  the  second  integral  vanishes  identically.  Therefore,  by  defining 


^s  “ 


I, 


^s+1 

^s 


I, 


^s+1 

Xs 


m(|)|  d| 
m(|)  d| 


(9) 


and 


ms=  j'''®-'^m(|)d|  (10) 

equation  (8)  becomes 

V 

cp{x)  =  cp^Xf)  +  J 
S=1 

where  cp^  = 

The  vanishing  of  the  second  integral  in  equation  (8)  when  |g  is  chosen  to  be  the 
interval  center  of  gravity  (c.g.)  suggests  that  this  is  an  optimum  location  of  the  point 
mass  to  minimize  the  error  bounds  (i.e.,  the  maximum  possible  error  involved  in 
replacing  the  continuous  system  by-the  lumped  mass  system).  Further  discussion  of 
these  approximation  errors  is  presented  in  a  subsequent  section. 

For  convenience,  xj  =  1^  is  chosen  to  be  the  reference  point  for  the  influence 
function  G(x,|).  Equation  (11)  for  x  =  |^,  where  r  =  1  to  v,  yields  a  set  of  equations 
which  can  be  conveniently  written  in  matrix  notation  as 


+  w 


(12) 
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where  Gj.  g  =  r,s  =  1  to  v. 

It  is  more  convenient  and  equally  correct  to  write  equation  (12)  as 

+  0)2  Gj.  g  nir  |^r^ 

where  r,s  =  1  to  v. 

Derivation  of  the  Lumped  Mass  Equations  for  a  Multibranch  Beam 

The  lumped  mass  approach  established  for  a  single  beam  is  extended  herein  to 
derive  similar  equations  of  motion  for  a  class  of  multibranch  beams;  this  class  consists 
of  a  main  beam  with  subsidiary  beams  attached  along  its  axis  like  a  tree  trunk  and  its 
branches.  Further,  the  subsidiary  or  branch  beams  are  parallel  to  the  main  branch. 

For  branch  identification  purposes,  the  main  branch  is  denoted  as  the  0th  branch; 
all  subbranches  are  denoted  by  the  mass  on  the  0th  branch  to  which  they  are  attached. 

The  lumped  masses  of  each  branch  are  numbered  consecutively  from  1.  The 
masses  of  the  various  branches  are  then  distinguished  from  each  other  by  a  branch  index 
which  appears  in  parentheses  after  the  mass  number.  The  same  mass-number — branch- 
number  index  system  is  also  used  for  the  x-locations  of  the  respective  masses. 

Consider,  for  a  moment,  the  main  branch  of  a  general  multibranch  beam.  This 
branch  can  be  divided  into  ''v(0)"  panels  or  intervals  of  arbitrary  length  as  shown  by  the 
dashed  line  in  sketch  3. 


Sketch  3 


The  total  weight  of  each  panel  can  then  be  concentrated  in  the  form  of  a  lumped 
mass  at  its  center  of  gravity  where  r  =  1  to  v(0).  (See  sketch  4.)  To  be  con¬ 

sistent  with  the  notation  introduced,  it  is  necessary  that  each  subbranch  attach  to  a  mass 
on  the  main  branch. 
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Sketch  4 


Then  as  shown  in  sketch  5,  any  subbranch  can  be  lumped  and  attached  to  the  appropriate 
mass  on  the  main  branch. 


Sketch  5 


Let  t5r(j),s(k) 

force  at  x  = 
as 


define  the  deflection,  relative  to  at  x  =  due  to  a  unit 

The  total  deflection  at  station  is  then  given  by  superposition 


v(0)  v(k) 

^r(3)=^l(0)+  2  I  ^r(j),s(k)  ^s(k) 
k=0  s=l 


(14) 


where 


is  the  rigid  body  translation  of  the  reference  point  x 


and  F 


IS 


^l(O)  lAgiu  uuuy  tJL  ctiiaicLLXuii  ui  tiic  1  cxcx  pvjiiix  cmvx  ^ 

the  axial  inertia  force  at  x  =  ^g(k)*  right-hand  summation  of  equation  (14)  is  over 

the  total  number  of  masses  of  a  given  branch  and  the  left-hand  summation  is  over  the 
total  number  of  branches. 


Equation  (14)  is  formulated  as  though  a  branch  is  located  at  every  mass  on  the 
0th  branch.  The  equation  holds  for  fewer  branches,  however,  if  the  subscripts  j  and  k 
only  assume  the  values  corresponding  to  branch  locations  (in  addition  to  0  for  the  main 
branch).  For  example,  if  branches  were  located  at  ^5(0)’  ^7(0)’ 


j,k  =  0,  1,  5,  7 
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The  subscripts  j  and  k  always  begin  with  0,  and  for  systems  consisting  only  of  the 
main  branch  they  only  assume  the  value  of  0. 

If  j  =  0  to  v(0)  and  equations  (3)  and  (4)  are  utilized,  equation  (14)  yields  a  set  of 
equations  which  can  be  written  in  the  following  matrix  form: 


+  O) 


^r(i),s(k) 


(15) 


where 


r  =  1  to  v(j) 
s  =  1  to  v(k) 
i,k  =  0  to  v(0) 

r  '\ 


‘^r(0) 

<Pr(l) 

l‘^r(2)f 


^r(j),s(k^ 


V.  •  J 

k - .- 


— 

^r(0),s(0) 

^r(0),s(l) 

®r(0),s(2) 

• 

• 

• 

‘^r(l),s(0) 

^r(l),s(l) 

^r(l),s(2) 

_ _  J 

• 

^r(2),s(0) 

^r(2),s(l) 

^r(2),s(2) 

L_  _  J 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 
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k 


rm 


"'s(O)^ 

0 

0 

u  _  _  _ _ _  _ 

i.  _  _ _ _  „ 

L 

0 

""3(1) 

0 

0 

O 

L 

''"'s(2)^ 

— 

• 

L 

r - 

• 

L 

• 

0 

• 

^  mm  imt  mm  mm  mm 

• 

r . 

u _ — 1  — 

- ^  1 

N 

— 

• 

^  _ _ 

< - 

Each  of  the  dotted-line  compartments  represents  a  partitioned  submatrix.  As  for  equa¬ 
tion  (12),  it  is  correct  and  useful  to  write  equation  (15)  as 


(16) 


Equations  (14),  (15),  and  (16)  are  then  respectively  analogous  to  equations  (11),  (12),  and 
(13)  which  were  previously  derived. 


Since  the  beam  is  considered  to  be  free  or  unconstrained,  longitudinal  equilibrium 
requires  that 


v(0)  v(j) 

I  I 


i=0  r=l 


=  0 


(17) 


Again,  by  utilizing  equations  (3)  and  (4)  equation  (17)  leads  to 

v(0)  v(3) 

I  I  ’’rO)  “r(j)  =  » 
j=0  r=l 

or  in  matrix  form  to 


(18) 


Premultiplying  equation  (16)  by 


1 

m  /.V 

r(]) 

,  utilizing  equation  (18),  and  solving  for 


^l(O) 
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1(0)  “irl*  “rO),  °r(j),s(k)||‘“r(i)jr''(» 


where 


”=  'IMf* 

_ )L  — 

is  the  total  mass  of  the  system.  Substituting  equation  (19)  into  equation  (16)  gives 


(j),s(k)  rr(i),s(k)  "^r(i)  If  r(j)J  ==  fr(i) 


where 


®r(j),s(k) 


'  H  -  D]|"r(j)/H 


The  index  notation  of  equation  (20)  has  utility  only  in  the  derivation  and  interpretation  of 
the  numbers  which  form  the  matrix  arrays.  Therefore,  to  simplify  the  notation  the  fol¬ 
lowing  matrices  are  redefined  and  used  interchangeably; 


®  ®a,b  ®r(j),s(k) 


G  =  G  b  - 


m  =  nio  =  m 


a  -  r“r(i) 


^  f  r(j) 


where 

a,b  =  1  to  n 
r  =  1  to  v(i) 
s  =  1  to  v(k) 
j,k  =  0  to  v(0) 

and  the  subscript  n  represents  the  total  number  of  lumped  masses  comprising  the 
system. 
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Equation  (20)  can  then  be  written  more  compactly  as 


where 


and 


(21) 


a,b  =  1  to  n 


— 

'v  - 

‘^a,b 

ma 

_  — 

X  = 


Equation  (21)  is  the  familiar  form  for  a  matrix  eigenvalue  problem  and  can  be 
solved  by  a  number  of  classical  procedures.  One  of  these  procedures  (an  iterative  solu¬ 
tion  involving  orthogonality  sweeping  for  determining  modes  higher  than  the  fundamental) 
is  presented  subsequently. 


Calculation  of  Influence  Coefficients 

Consider  a  beam  fixed  at  station  with  a  unit  load  applied  in  the  positive 

x-direction  at  station  as  shown  in  sketch  6. 


Sketch  6 

The  deflection  at  any  point  x^q^  is  then  defined  as  an  influence  coefficient.  These 
deflections  or  influence  coefficients  are  derived  for  all  combinations  of 

j  2  k  through  the  use  of  elementary  one-dimensional  beam  theory,  as  follows: 
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For  j  =  k  =  0, 


For  j  =  0  and  k  #  0, 


For  j,k  it  0, 


,s(0)  j! 


r(0)  dx 


X  ,  ,  A(0)  E(0) 

Xi(o) 


(r  5  s)  (22) 


^r(0),s(k)  =, 

y^r(0) 

^""1(0) 

dx 

A(0)  E(0) 

(r 

£k) 

(23) 

®r(0),s(k)  =, 

J^k(0) 

dx 

A(0)  E(0) 

(r 

>k) 

(24) 

^1(0) 

dx 

(i 

<k) 

(25) 

®r(i),s(k) 

Jx 

^l(O) 

A(0)  E(0) 

G  =  r^j(O)  dx  r^r(j)  dx  /26) 

^r(j),s(k)  -  A(0)  E(0)  A(j)  E(j) 

(j  =  k  and  r  ^  s) 

The  absolute  value  of  the  integral  in  equation  (26)  allows  for  the  possibility  of  a  branch  in 
either  the  positive  or  negative  x-direction. 

All  values  of  the  influence  coefficients  for  j  >  k  can  be  determined  from  Maxwell's 
reciprocal  theorem;  thus, 

®r(i),s(k)  ""  S(k),r(j) 

The  1/AE  function  in  the  Gj.(j)^s(l5;)  integral  is  generally  a  highly  discontinuous 
function  and,  therefore,  a  numerical  form  of  integration  is  expedient.  A  simple  quadra¬ 
ture  numerical  integration  will  give  excellent  results  and  yields  the  exact  answer  provided 
that  the  1/AE  curve  is  a  straight  line  between  intermediate  integration  points.  For  a 
curved  line,  the  intermediate  integration  points  should  be  generously  chosen  so  that  the 
curve  can  be  approximated  as  a  straight  line  between  them. 


Application  of  the  Method  to  Torsional  Vibrations 

The  Saint-Venant  torsion  theory  for  prismatical  bars  leads  to  the  following  classical 
equation  of  motion  for  torsional  vibrations: 
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(28) 


where 
C  =  C(x) 
I  =  I(x) 

9  =  e{x,t) 


a 


torsional  rigidity  (ref.  9) 

polar  mass  moment  of  inertia  per  unit  length 

cross-section  rotation 


For  a  circular  cross  section  the  torsional  rigidity  is 

C  =  JG  (29) 

where  J  is  the  polar  moment  of  inertia  of  the  cross  section  and  G  is  the  shear 
modulus.  Equation  (28),  for  bodies  of  revolution,  thus  becomes 


The  similar  one-dimensional  equation  of  motion  for  longitudinal  vibrations  on  which 
the  derivations  of  foregoing  sections  are  founded  is 


where 


A  =  A(x) 
E 

m  =  m(x) 
u  =  u(x,t) 


cross-sectional  area 
modulus  of  elasticity 
mass  per  unit  length 
cross-section  displacement 


If  equation  (28)  (or  eq.  (30))  is  compared  with  equation  (31),  it  is  seen  that  C 
(or  JG)  is  analogous  to  AE,  I  is  analogous  to  m,  and  6  is  analogous  to  u.  There¬ 
fore,  by  replacing  AE,  m,  and  u  with  C,  I,  and  6,  respectively,  all  the  previously 
derived  equations  for  the  longitudinal  vibration  of  beams  can  be  applied  directly  to  tor¬ 
sional  analyses  of  structures  whose  equation  of  motion  is  compatible  with  equation  (28). 
This  compatibility  is  only  satisfied  when  the  elastic  axes  of  all  branches  are  everywhere 
concentric. 
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Derivation  of  Orthogonality  Relationships 

Several  of  the  classical  iteration  techniques  for  solving  equation  (21)  depend  upon 
orthogonality  relationships  for  obtaining  modes  higher  than  the  fundamental.  These  rela¬ 
tionships  are  derived  herein. 

From  equation  (16),  the  equation  for  the  pth  natural  mode  can  be  written  as 


jjp(p) 

where  <pj(p)  =  <Pj(q)(p)-  Premultiplying  equation  (32)  by  |^(q^|^n^,  where  [5>(q^  is 
the  qth  natural  mode,  gives 

[^(q^[m^{^3(p)}  =  [5p(q)J  ^{l}  (pj(p)  +  w^(p)[f  (q^|jn^  @1^  ^(p)}  (33) 

Proceeding  in  a  similar  manner  for  the  qth  mode  gives 

l^(p)J  If 


^i(p)fi  1+  ^^(p) 


rm 

m 

- 


(32) 


The  transposition  of  equation  (18)  gives 

If  (q)J  ^  if  ^pU  KI(^}  =  ® 


Also,  since 


[G^  and  [in^  are  symmetrical, 


Lf(q)J|ni][G]Im]{<p(p)}  =  [5)(p)J[n^[(:^[m]{^(q)) 


(36a) 


and 


l^p(q)|  [m]  •(f  (P)}  =  If  (pUl^m]  <(f  (q)}  (36b) 

Therefore,  subtracting  equation  (34)  from  equation  (33)  and  utilizing  equations  (35)  and 
(36)  gives 

(w^(p)  -  aj^(q)]  |^<?(q)j  [m]  [g]  [m]  <{^(p)}  =  0  (37) 

Then  from  equation  (37)  for  p  #  q 

If  (q)J  [m]  [pj  ^{^(P)}  =  0  (36a) 
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Substituting  equations  (35)  and  (38a)  into  equation  (33)  gives 

(q^  [m^  (p)}  =  0  (38b) 

for  p  q. 

Equations  (38)  are  then  two  valid  orthogonality  relationships.  Equation  (38b)  is 
generally  preferred,  however,  because  of  its  simplicity. 

Solutions  to  Eigenvalue  Problem 

It  is  beyond  the  scope  of  this  paper  to  enumerate  the  iterative  procedures  for 
solving  eigenvalue  problems  and  to  discuss  the  advantages  of  each.  The  methods  can 
generally  be  categorized  into  two  basic  types:  those  which  replace  a  trial  vector  with  an 
improved  trial  vector,  and  those  which  replace  a  square  matrix  with  an  improved  square 
matrix.  A  method  of  each  type  is  discussed  briefly  herein. 

Iteration  and  sweeping.-  A  classical  method  for  solving  eigenvalue  problems 
involves  the  replacement  of  a  trial  vector  by  an  improved  vector  obtained  by  multiplying 
the  trial  vector  by  the  coefficient  matrix.  The  method  is  classical  and  is  documented  in 
numerous  works  (e.g.,  ref.  10).  Applying  this  process  to  equation  (21),  however,  gives 
only  the  fundamental  mode. 

Equation  (21)  can  be  modified  by  any  one  of  several  different  procedures  to  deter¬ 
mine  higher  modes.  One  particular  method  which  has  provbd  very  satisfactory  is  out¬ 
lined  in  reference  11.  By  utilizing  the  orthogonality  relationship  (eq.  (38b)),  the  higher 
modes  can  be  determined  by  the  iteration  of  the  following  modified  matrix  equation: 

|P(P)]{‘^}  =  ^  {^}  (P  =  2,  3,  .  .  .,n-l)  (39) 


where 


[p(p)]  =  [D][s(pj] 

is  defined  by  equation  (21) 


ma 
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and 


ii(p)] = 


Cs(P  - 1)]  - 


(pjjp  - 1) 

ST  “ 

ma 

(pji'p  - 1) 

ma 

L 

-«1 

(p  =  3,  4,  5,  .  .  n) 


The  [S(p)J  matrix  is  a  so-called  sweeping  matrix  which  removes  the  components 
of  all  modes  below  the  pth  mode  from  the  trial  iteration  vector  This  method  is 

suitable  for  desk  calculator  or  digital  computer  solution. 

The  characteristic  polynomial  of  equation  (39)  clearly  is  of  order  "n";  therefore, 
the  system  has  "n"  eigenvalues  (X's).  Only  n  -  1  of  these  eigenvalues  and  modes  can 
be  obtained  by  iteration.  The  remaining  mode  is  spurious  and  corresponds  to  a  zero 
eigenvalue  (X).  The  existence  of  this  spurious  mode  can  be  proved  as  follows: 

It  is  noted  in  equation  (21)  that 

[d]  =  [B][G][ni]  (40) 

Since  [g]  is  the  influence  coefficient  matrix  for  a  beam  fixed  at  the  first  mass,  the  first 
row  and  first  column  of  [g]  are  null,  and  (g^  is,  therefore,  singular.  It  then  follows 
from  the  theory  of  determinants  (ref.  12,  p.  25)  that  [j^  is  also  singular  —  that  is. 


[B] 


=  0 


(41) 


It  can  be  further  concluded  from  the  theory  of  determinants  (ref.  12,  pp.  68-69)  that  the 
singularity  of  [d]  necessitates  a  zero  eigenvalue  (X)  of  equation  (21). 

The  mode  corresponding  to  this  zero  eigenvalue  must  satisfy  equation  (21)  for 
X  =  0.  Thus, 


=  0  (42) 

The  singularity  of  [d^  is  a  necessary  and  sufficient  condition  to  guarantee  a  nontrivial 
solution  to  equation  (42)  (ref.  13).  This  nontrivial  solution  has  no  physical  significance 
and  can  be  shown  by  substitution  to  be 
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(43) 


0 


From  physical  considerations  there  are  n  -  1  elastic  modes  and  one  rigid  body 
(zero  frequency)  mode  associated  with  the  lumped  mass  system  which  led  to  equation  (21). 
The  transformation  resulting  from  equation  (19)  eliminated  the  rigid  body  mode  and 
replaced  it  by  the  aforementioned  spurious  mode.  Iteration  of  equation  (21)  thus  only 
yields  the  n  -  1  dominant  elastic  modes. 

Matrix  diagonalization.-  The  computer  age  has  revived  interest  in  methods  for 
solving  characteristic  value  problems  in  which  the  coefficient  matrix  is  brought  to  diag¬ 
onal  or  nearly  diagonal  form  by  a  sequence  of  orthogonal  transformations.  Among  these 
are  many  variations  of  an  old  procedure  formulated  by  Jacobi  in  which  the  matrix  is  diag¬ 
onalized;  this  matrix  diagonalization  leads  to  solution  for  all  modes  and  frequencies 
simultaneously.  In  a  recent  innovation  by  Givens  (ref.  10,  pp.  336-340)  the  original 
matrix  is  transformed  to  tridiagonal  form  and  the  eigenvalues  are  then  computed  from 
a  Sturm  sequence. 

All  such  methods,  however,  require  the  eigenvalue  problem  to  be  in  the  following 

form: 


[£](?,}  =  (44) 

where  [e^  is  a  real  symmetric  matrix  and  A  is  the  eigenvalue.  Equation  (21)  can  be 
put  in  the  form  of  equation  (44)  by  the  following  procedure. 

A  sweeping  matrix  similar  to  that  of  equation  (39)  is  utilized  to  remove  the  rigid 
body  mode  -(^(0)^  =  ^1^  from  the  trial  vector.  This  sweeping  matrix  is 


[S(1)]  =  H  - 


(45) 


where  is  the  transpose  of  the  [b^  matrix  of  equation  (20).  Since  the  rigid  body 

mode  has  already  been  eliminated,  this  step  is  redundant.  It  does,  however,  lead  to  the 
desired  symmetrical  form. 

It  can  be  seen  from  equation  (18)  that  [S(l)]  possesses  the  property 

[S(l3{<p(p)}  =  fyip)}  (p  =  1,  2,  3,  .  .  .,  n  -  1)  (46) 
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(47) 


Substituting  equations  (45)  and  (46)  into  equation  (21)  gives 


Then,  introducing  the  change  in  variables 


into  equation 


Equation  (49) 


(y)  = 

(48) 

(47)  leads  to 

(x) 

(49) 

is  then  in  the  form  of  equation  (44)  where 

CE]=Ml/2CB]CGW[m]l/2 

After  solving  equation  (49)  for  {y} ,  the  original  modal  matrix  is  obtained  from  equa¬ 
tion  (48)  as 

(50) 

The  system  given  by  equation  (47)  also  has  a  spurious  mode  which  corresponds  to  a 
zero  eigenvalue  (X).  The  proof  is  similar  to  the  one  given  in  the  preceding  section.  The 
spurious  modal  vector  in  this  case,  however,  turns  out  to  be  identical  to  the  rigid  body 
mode  in  contrast  to  the  form  of  equation  (43).  This  identity  can  be  shown  as  follows: 

From  equation  (45),  equation  (47)  can  be  written  as 

where  X  =  If  X  =  0,  the  spurious  mode  must  satisfy 

Cb][g|M(s(i)]  {(py = 0 

This  equation  is  obviously  satisfied  by  the  rigid  body  mode  =  {l}’  since 


[S(l3{l}  =  0 
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NUMERICAL  EXAMPLES 


Idealized  Beam 

For  the  purpose  of  showing  the  operations  of  the  matrices  derived  in  the  foregoing 
section,  the  idealized  beam  shown  in  sketch  7  is  considered.  This  example  consists  of  a 
solid  uniform  beam  attached  to  the  inside  of  a  uniform  cylindrical  shell  by  a  massless, 
infinitely  rigid  point  attachment  which  coincides  with  the  shell  center  of  gravity. 


m,AE 


n:' . ^ . ^  . V . V  V  \  V  V  V  ^  v  ^ 

^^*.^2m,2AE 

////// 

/ / / / 

.  L  J 

r  2 

3  1 

Sketch  7 

For  simplicity  of  illustration,  the  cylindrical  shell  is  represented  by  three  equal 
panels  and  the  solid  inner  beam  by  one  panel.  The  total  mass  of  each  panel  can  be  con¬ 
sidered  concentrated  at  the  panel  center  of  gravity  as  shown  by  the  analogous  beam  in 
sketch  8. 


2L/3 


m 


1(0) 


m 


5L/6 


2(0) 


m 


3(0) 


Sketch  8 

Assume  that  the  cross-sectional  properties  of  the  cylindrical  shell  are  m(x)  and 
A(x)  E(x)  and  those  of  the  solid  inner  beam  are  2m(x)  and  2A(x)  E(x).  The  matrices 
used  to  formulate  equation  (21)  are  then 
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1  0  0 
M  0  0  1 


0  0  0 


[G]  = 


0  0  0  10 
0  4  4  I  4 


12AE  0  4  8 


0  4  4  1  5 

I 


10  0  0 
0  10  0 
0  0  10 
0  0  0  1 


1111  1000 
1  1  1  1  1  0  10  0 
5  1111  0  0  10 
1  1  1  1  0  0  0  2 


4  -1  -1  -2 
1-1  4-1  -2 

5-1-1  4  -2 

-1  -1  -1  3 


where  the  0th  branch  parameters  are  indicated  by  the  upper  left-hand  partition  matrices. 
Substituting  equations  (51),  (52),  and  (53)  into  equation  (21)  gives 


0  -8  -10  -18 
niL^O  2  0  2  = 

90AE  0  2  10  2  I  (  “ 

0  2  0  7  I 


After  iteration  of  equation  (54)  the  following  solutions  are  obtained  for  the  funda- 
mental  flexible  mode: 


-a(«  =  I 


Substituting  equations  (55)  into  the  expression  for  [^(2)]  of  equation  (39)  gives 
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-N 

-1 

1-1  0  1  qj 

10  0  0 

0 

0  10  0 

11 

0  0  10 

UJ 

0  0  0  2 

0 

0  1 

0 

1-1  0  1  qj 

1 

0 

0 

0 

r-i] 

0 

0  0 

1 

0 

1 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

2 

0 

L  J 

10  10 
1  0  2  0  0 
2  10  10 
0  0  0  2 


(56) 


Substitution  of  equation  (56)  into  the  expression  for 


[p(p^  of  equation  (39)  gives 


0  -8  -10  -18 

10  10 

mL^ 

0  2  0  2 

0  2  0  0 

180AE 

0  2  10  2 

10  10 

0  2  0  7 

0  0  0  2 

mL^ 

90AE 


-5  -8  -5 

0  2  0 

5  2  5 

0  2  0 


-18 

2 

2 

7 


Substituting  equation  (57)  into  equation  (39)  gives 


mL^ 

90AE 


-5  -8  -5  -18 

C 

ri 

0  2  0  2 

1 

^  1 

5  2  5  2 

< 

“^a 

11 

^a 

0  2  0  7 

L  J 

L  J 

(57) 


(58) 


After  iteration  of  equation  (58),  the  following  solutions  are  obtained  for  the  second  flex¬ 
ible  mode: 
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r 

k(2)i 

V  J 


1.0000 

0.2986 

1.0000 

0.8507 


0,(2)  = 


-(59a) 


(59b) 


In  a  similar  manner,  the  solutions  for  the  third  flexible  mode  are  determined  to  be 


(^0.1489 
1.0000 
-0.1489 
-0.351^ 


(60a) 


0,(3) 


8.33,51 

L  riir 


(60b) 


A  plot  of  the  three  flexible  modes  obtained  in  this  example  are  shown  in  figure  1. 

It  is  interesting  to  note  that  even  with  the  inclusion  of  branches,  the  number  of  nodal 
points  corresponds  to  the  flexible  mode  number  (i.e.,  the  first  flexible  mode  has  one  node 
point,  etc.). 


Typical  Launch  Vehicle 

To  show  the  application  of  the  matrix  solution  to  a  more  realistic  beam,  a  numerical 
example  of  an  application  to  a  solid-fuel  booster  system  is  presented.  The  vehicle  param¬ 
eters  are  typical  of  those  of  the  Scout  class.  The  payload  parameters  were  selected  for 
purposes  of  illustration  and  are  not  necessarily  representative  of  an  existing  payload. 

The  following  assumptions  are  employed: 

(1)  The  solid-fuel  mass  is  considered  to  adhere  to  the  sides  of  the  vehicle  along  the 
length  and  has  no  motion  relative  to  the  vehicle. 

(2)  The  contribution  of  the  fuel  stiffness  to  the  vehicle  axial  stiffness  is  negligible. 

(3)  Damping  is  considered  to  be  negligible. 

(4)  All  deformations  are  one-dimensional  and  no  consideration  is  given  to  bending 
or  breathing  effects  of  the  cylindrical  shell  walls. 

Figures  2  and  3  show  the  rocket-vehicle  mass  and  axial-extensional-coefficient 
data,  respectively,  in  a  desirable  form  for  the  analyst.  The  data  for  these  graphs  are 
given  in  tables  1,  2,  and  3.  The  analogous  discrete  mass  system  is  shown  in  figure  4. 

The  X  origin  was  chosen  to  coincide  with  the  first  discrete  mass  rather  than  the  end  of 
the  beam. 
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The  first  three  normalized  natural  modes  and  their  corresponding  natural  frequen¬ 
cies  obtained  from  the  iteration  of  equation  (21)  are  shown  in  figure  5.  Note  the  high 
amplitude  of  the  payload  motion  in  both  the  first  and  second  modes  of  vibration  as  a 
result  of  the  payload  frequency  being  close  to  the  lower  frequencies  of  the  vehicle.  Such 
a  design  can  cause  extensive  acceleration  loading  on  sensitive  payload  instrumentation 
and  should  be  avoided  if  at  all  possible. 

In  both  the  first  and  second  modes,  the  vehicle  is  deflected  in  a  manner  analogous  to 
that  of  the  fundamental  vehicle  mode  without  the  payload.  The  addition  of  the  payload 
branch  introduces  a  whole  range  of  secondary  coupled  frequencies  in  which  the  payload 
oscillations  are  either  in  phase  or  out  of  phase  with  the  vehicle  oscillations. 

DISCUSSION  OF  COMPUTATIONAL  ACCURACY 


The  previously  derived  equations  pertinent  to  the  solution  were  programed  for  the 
IBM  7094  electronic  data  processing  system.  Some  of  the  points  of  interest  in  the  appli¬ 
cation  of  this  method  are  discussed. 

Choice  of  the  Number  of  Lumped  Masses 

The  matrix  solution  used  herein  is  subject  to  an  inherent  error  dependent  upon  the 
validity  of  an  analogous  lumped  mass  approximation  to  a  continuous  beam.  To  gain  an 
insight  into  the  accuracy  of  the  matrix  solution,  comparisons  were  made  between  the 
theoretical  and  approximate  results.  As  mentioned  previously,  closed  form  solutions  for 
beams  of  arbitrarily  varying  cross  sections  are  very  difficult  to  obtain.  Therefore,  com¬ 
parisons  were  limited  to  beams  with  uniform  and  exponentially  varying  cross  sections. 
For  the  approximate  solution  the  beams  were  divided  into  n  equally  spaced  intervals 
with  the  lumped  mass  at  the  center  of  gravity.  The  exact  solutions  for  these  beams  are 
derived  in  appendix  A. 

A  standard  accuracy  parameter  for  frequency  comparison  is  the  percent  error 
which  is  calculated  as  follows: 

Percent  error  = - ^ - x  100 

^exact 

No  convenient  standard  exists  for  mode  comparison.  For  the  purposes  of  this 
paper,  a  root- mean- square  percent  error  (derived  in  appendix  B)  is  used. 

The  percent  error  of  the  first  five  elastic  modes  as  a  function  of  the  total  number 
of  lumped  masses  n  is  given  in  figures  6  and  7  for  the  frequencies  and  mode  shapes, 
respectively.  It  can  be  seen  in  both  figures  that  the  cross-sectional  exponential  factor  jS 
has  very  little  effect  on  the  percent  error,  although  the  range  of  /3  considered  produced 


a  wide  variation  in  the  distribution  of  cross-sectional  area.  These  results  indicate  that 
both  frequency  and  mode  accuracy  depend  more  on  the  number  of  masses  than  on  the 
variation  of  cross-sectional  area.  As  expected,  the  error  curve  asymptotically 
approaches  zero  for  high  values  of  n.  In  figure  7,  considerable  slope  discontinuities  are 
apparent  in  the  curves  toward  the  origin.  These  discontinuities,  for  a  given  mode,  occur 
in  the  rar^e  where  n  is  either  1  or  2  greater  than  the  mode  number.  It  is  noted  that 
in  this  range  of  n  the  addition  of  an  extra  mass  may  or  may  not  be  detrimental  to  the 
mode  shape. 

For  most  engineering  applications,  a  1-percent  error  limit  is  sufficient.  There¬ 
fore,  the  1-percent-error  data  for  the  first  five  modes  from  figures  6  and  7  are  pre¬ 
sented  in  figure  8.  It  can  be  seen  that  the  solid  line  corresponding  to  the  equation 

n  =  6i  +  1 


(where  i  is  the  highest  mode  to  be  determined)  gives  a  "rule  of  thumb"  for  determining 
an  adequate  number  of  masses  to  obtain  both  the  frequencies  and  the  modes  with  this 
accuracy. 


Advantages  of  Influence  Coefficient  Approach 

In  discrete  element  structural  analysis,  two  approaches  can  be  taken:  the  stiffness 
(spring)  approach  and  the  flexibility  (influence  coefficient)  approach.  The  tendency  is  to 
adopt  the  stiffness  approach  without  consideration  of  accuracy  loss  because  the  stiffness 
matrix  is  generally  simpler  to  formulate  than  the  flexibility  matrix.  The  stiffness  matrix 
approach,  however,  necessitates  the  formulation  of  the  matrix  eigenvalue  problem  such 
that  the  eigenvalue  corresponding  to  the  highest  natural  frequency  is  dominant. 

The  prime  advantage  of  the  flexibility  approach  is  that  the  problem  can  be  formu¬ 
lated  so  that  the  mode  corresponding  to  the  lowest  natural  frequency  is  dominant.  It  is 
generally  accepted  that  the  various  means  of  solving  eigenvalue  problems  lose  accuracy 
when  proceeding  from  the  dominant  mode  to  the  subdominant  modes.  The  magnitude  of 
this  accuracy  degeneration  varies  somewhat  with  different  methods  of  solution,  but  the 
author's  experience  and  that  of  Batchelder  and  Wada  (ref.  4)  indicate  that  the  error  is 
dependent  on  the  ratio  of  the  dominant  eigenvalue  to  the  subdominant  eigenvalue.  Also, 
experience  has  shown  that  as  this  ratio  becomes  large,  accuracy  is  lost  in  the  subdomi¬ 
nant  modes. 

A  sizable  spread  between  the  dominant  and  subdominant  frequencies  can  result  from 
a  large  number  of  degrees  of  freedom  and/or  large  variations  in  the  stiffness  or  mass  in 
different  parts  of  the  system.  The  stiffness  matrix  approach  can  lead  to  complete  loss  of 
accuracy  in  the  lower  modes,  which  are  of  prime  interest;  whereas,  the  flexibility  matrix 
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approach  only  loses  accuracy  in  the  less  important  higher  modes.  This  fact  is  simply 
illustrated  by  comparative  solutions  of  the  spring-mass  system  shown  in  sketch  9. 


Sketch  9 

From  sketch  9  it  can  be  seen  that  as  the  value  of  77  becomes  large,  k4  becomes 
large  and  the  lower  frequencies  of  the  system  approach  those  of  the  system  shown  in 
sketch  10(a).  In  addition,  the  value  of  k3  becomes  negligible  in  comparison  with  the 
value  of  k4,  and  the  highest  frequency  of  the  system  approaches  the  frequency  of  the 
uncoupled  two  mass  system  shown  in  sketch  10(b).  The  frequencies  of  these  two  simple 
degenerate  systems  can  be  determined  with  high  accuracy  (by  any  method)  and  are, 
therefore,  used  as  a  reference  for  a  subsequent  comparative  study. 


(b)  Fourth  mode. 

Sketch  10 

The  natural  frequencies  were  obtained  as  a  function  of  increasing  values  of  r]  for 
the  system  of  sketch  9,  by  methods  employing  both  the  flexibility  and  stiffness  matrices. 
The  ratios  of  these  frequencies  to  the  aforementioned  reference  frequencies  are  given  in 
table  4.  If  there  were  no  loss  of  accuracy  in  computation,  the  ratios  would 

approach  unity  as  77  becomes  large.  A  departure  from  this  pattern  indicates  loss  of 
numerical  accuracy  in  computation. 

It  can  be  seen  that  method  I  (Givens  solution  to  stiffness  matrix)  gives  good  answers 
up  to  the  point  77  =  6,  which  corresponds  to  an  approximate  ratio  of  10^  between  the  high¬ 
est  and  the  lowest  mode  frequencies.  Beyond  this  point  there  is  complete  accuracy  loss 
in  all  but  the  highest  mode.  Method  II  (iteration  and  sweeping  solution  to  stiffness  matrix) 
indicates  similar  accuracy  loss  in  the  lower  modes;  although,  the  drop-off  is  not  as  sharp 
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as  in  method  I.  Method  III  (iteration  and  sweeping  solution  to  flexibility  matrix)  and 
method  IV  (Givens  solution  to  flexibility  matrix)  give  very  similar  answers  and  show  no 
accuracy  loss  except  in  the  highest  mode.  For  low  values  of  ry  all  methods  are  seen 
to  give  nearly  the  same  answers. 

This  example  does  not  represent  a  real  structure  but  it  does  point  out  the  type  of 
problems  that  can  occur  and  have  occurred  to  the  author's  knowledge  in  practical  analy¬ 
ses  of  structures  with  widely  varied  mass  and  elastic  properties.  It  also  clearly  demon¬ 
strates  the  increased  lower  mode  computational  accuracy  obtained  by  the  flexibility 
approach. 

In  many  structural  problems,  the  spread  between  the  highest  and  lowest  natural 
frequencies  is  small  enough  to  yield  satisfactory  results  by  the  stiffness  matrix  approach. 
There  is  always  the  element  of  doubt,  however,  and  the  analyst  must  be  constantly  aware 
of  impending  accuracy  loss  in  the  lower  modes.  The  influence  coefficient  approach  of 
this  paper  is  thus  advocated  to  yield  confidence  in  the  lower  modes  and  to  alleviate  the 
extra  work  associated  with  an  approach  which  at  times  breaks  down. 

CONCLUDING  REMARKS 

An  influence  coefficient  matrix  approach  to  the  problem  of  longitudinal  unconstrained 
vibrations  of  nonuniform  multibranch  beams  is  presented.  The  derivation  is  completely 
general  and  requires  only  basic  parameter  changes  from  one  application  to  another.  By 
making  a  change  in  variables  the  method  can  also  be  applied  to  torsional  vibration  prob¬ 
lems.  The  derived  matrix  equations  have  been  programed  for  the  IBM  7094  electronic 
data  processing  system  and  have  given  excellent  results  in  numerous  applications. 

A  complete  numerical  example  is  given  to  illustrate  the  computations  involved  and 
a  detailed  illustration  of  the  method  applied  to  a  launch  vehicle  is  provided. 

Comparison  of  the  method  with  available  exact  solutions  indicates  that  the  accuracy 
of  the  solution  is  practically  independent  of  cross-sectional  variations  but  is  contingent 
upon  the  chosen  number  of  lumped  masses  in  the  analogous  system.  An  approximate  rule 
is  given  for  determining  the  number  of  lumped  masses  to  yield  less  than  a  1-percent 
error  in  frequencies  and  mode  shapes. 

The  influence  coefficient  matrix  method  is  also  compared  with  the  stiffness  matrix 
approach  and  is  shown  to  be  more  accurate  in  the  lower  modes  when  there  is  a  wide  var¬ 
iation  in  the  magnitudes  of  the  frequencies  of  the  system. 

Langley  Research  Center, 

National  Aeronautics  and  Space  Administration, 

Langley  Station,  Hampton,  Va.,  September  19,  1966, 

124-11-05-30-23. 
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APPENDIX  A 


CLASSICAL  WAVE  EQUATION  SOLUTION  FOR  BEAM  WITH 
EXPONENTIALLY  VARYING  CROSS  SECTION 


For  completeness,  the  solution  for  a  beam  with  exponentially  varying  cross  section 
is  given.  The  classical  wave  equation  for  longitudinal  vibrations  is 


—  (AE  — 
axV  9x 


A 

=  pA-^ 

dr 


where 


A  =  A(x)  cross-sectional  area 

E  modulus  of  elasticity 

p  =  p(x)  mass  density 

u  =  u(x,t)  displacement 

For  undamped  simple  harmonic  motion  of  frequency  o), 

u(x,t)  =  cp{x)e^^^ 

Substituting  this  expression  into  equation  (Al)  gives 

The  beam  under  consideration  has  a  cross-sectional  area  such  that 

A(x)  =  A(0)e^^/^ 

where  L  is  the  beam  length  and  (3  is  a  constant. 

Substituting  this  expression  into  equation  (A2)  gives 

=  0  (A3) 

2 

j3  ,  the  solution  to  equation  (A3)  is 


30 


APPENDIX  A 


<p(x)  =  cos 

where  A  and  B  are  constants  of  integration. 

If  the  free-free  boundary  conditions 

^(0)  =  ^(L)  =  0 

dx  dx 

are  utilized  to  evaluate  A  and  X,  equation  (A4)  becomes 

,y(x)  =  Ce-'^‘-(cos  =  ^jlsin!s)  (A5) 

where 

C  =  an  arbitrary  constant 

n  =  0,  1,  2,  ,  .  .  an  integer 


^  +  B  sin  ^ 
L 


Normalizing  the  displacement  on  (p(0)  gives 


<?(x)  ^ 

<p(0) 


nTPc 

L  nTT 


sin 


n7rx\ 

L/ 


The  frequency  of  vibration  is  determined  from  the  evaluation  of  X  to  be 


(A6) 


(A7) 


Equations  (A6)  and  (A7)  give  the  theoretical  solution  for  the  natural  modes  and  frequen¬ 
cies  of  vibration.  These  equations  reduce  to  the  classical  solution  for  a  uniform  beam 
if  j3  =  0. 
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APPENDIX  B 


ROOT-MEAN- SQUARE  PERCENT  ERROR 


To  compare  an  approximate  and  theoretical  curve  by  a  single  meaningful  parameter 
is  difficult  since  each  curve  consists  of  an  infinite  number  of  points.  This  difficulty  is 
compounded  by  the  more  basic  problem  of  how  to  initially  aline  the  two  curves. 

The  approximate  mode  shapes  are  normalized  by  an  arbitrary  constant.  Therefore, 
adjusting  this  constant  such  that  the  error  between  the  theoretical  and  approximate  curves 
is  a  minimum  yields  a  logical  approach  to  alinement. 

One  method  for  determining  an  error  parameter  is  to  average  the  percent  error  of 
the  points  along  the  curve  or  to  take  the  maximum  value.  This  process,  however,  leads 
to  a  distorted  accuracy  picture  due  to  the  high  percent  error  in  the  nodal  point  regions. 

A  somewhat  abstract  but  yet  more  realistic  approach  is  to  base  the  percent  error  on  the 
average  value  of  the  theoretical  curve.  The  derivation  of  such  a  parameter  is  included 
herein. 

Let  E^  represent  the  mean  square  error  between  a  typical  theoretical  curve  and 
an  approximate  curve  as  depicted  in  sketch  11, 


X 


Sketch  11 


The  value  of  E^  is  then  given  as 


N?(x)  -  cp{x) 


dx 


(Bl) 


where  N  is  an  arbitrary  normalizing  coefficient  of  ?(x).  For  the  minimum  value 
of  e2 


d(E^) 

dN 


=  0 


(B2) 


Performing  the  differentiation  of  equation  (Bl)  indicated  by  equation  (B2)  and  solving  for 
N  gives 
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Equation  (B5)  has  been  programed  for  the  IBM  7094  electronic  data  processing  sys¬ 
tem.  Since  the  curve  ^(x)  was  defined  by  a  limited  number  of  points,  a  second-order 
interpolation  scheme  was  used  for  the  numerical  integrations  of  the  integrals  in  equa¬ 
tions  (Bl),  (B3),  and  (B4). 
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TABLE  1.-  VEHICLE  MASS  DISTRIBUTION 


(a)  U.S.  Customary  Units 


X, 

in. 

2 

Ibf-sec'^ 

in2 

in. 

Ibf-sec^ 

in^ 

X, 

in. 

P 

Ibl-sec*^ 

in^ 

-15.8 

0.038808 

458.9 

0.004788 

496.299 

0.051756 

0 

.046624 

468.449 

.004788 

497.299 

.051756 

19.1 

.046624 

468.449 

.027562 

497.299 

.045316 

19.1 

.071398 

468.949 

.027562 

520.789 

.045316 

29.6 

.071398 

468.949 

.019565 

520.789 

.046553 

29.6 

.049883 

469.199 

.019565 

522.199 

.046558 

102.2 

.049883 

469.199 

.021480 

522.199 

.004710 

102.2 

.059536 

469.389 

.021480 

526.789 

.004710 

173.3 

.059536 

469.389 

.023784 

526.789 

.0159318 

173.3 

.100380 

469.949 

.023551 

531.099 

.0159318 

177.3 

.100380 

469.949 

.028338 

531.099 

.027685 

177.3 

.007668 

470.949 

.028080 

533.049 

.027685 

186.4 

.007668 

470.949 

.034161 

533.049 

.005435 

186.4 

.018549 

471.199 

.034032 

537.269 

.00  54  3  5 

188.4 

.018549 

471.199 

.018401 

537.269 

.005823 

188.4 

.012378 

474.949 

.017728 

538.809 

.005564 

195.9 

.012378 

474.949 

.011232 

538.809 

.006185 

195.9 

.025331 

476.199 

.011077 

540.199 

.005900 

1  207.4 

.025331 

476.199 

.006315 

545.199 

.005176 

1  207.4 

.023834 

476.749 

.006211 

548.309 

.004790 

j  214.9 

.023834 

476.749 

.008463 

548.309 

.008851 

214.9 

.024559 

482.285 

.008230 

552.0 

.00880 

317.1 

.024559 

482.285 

.014880 

554.059 

.008644 

317.1 

.044689 

482.749 

.014880 

554.059 

.001941 

322.1 

.044689 

482.749 

.012681 

559.949 

.001579 

322.1 

.013808 

483.199 

.012681 

559.949 

.000924 

331.4 

.013808 

483.199 

.011439 

562.0 

.00070 

331.4 

.025466 

485.999 

.011439 

565.199 

.000655 

345.4 

.025466 

488.285 

.011827 

566.0 

.000500 

345.4 

.022530 

488.285 

.005952 

573.46 

.000466 

350.4 

.022530 

488.785 

.006315 

573.46 

.001633 

350.4 

.024559 

488.785 

.018245 

577.629 

.001633 

451.9 

.024559 

490.799 

.019358 

577.629 

.001568 

451.9 

.044689 

490.799 

.045031 

577.629 

0 

456.9 

.044689 

496.299 

.045031 
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TABLE  1.-  VEHICLE  MASS  DISTRIBUTION  -  Concluded 

(b)  SI  Units 


-0.4013 

0 

.4851 

.4851 

.7518 

.7518 

2.5959 

2.5959 

4.4018 

4.4018 

4.5034 

4.5034 

4.7346 

4.7346 

4.7854 

4.7854 

4.9759 

4.9759 

5.2680 

5.2680 

5.4585 

5.4585 
8.0543 
8.0543 
8.1813 
8.1813 
8.4176 
8.4176 
8.7732 
8.7732 
8.9002 
8.9002 


m, 

N-sec^ 

X, 

m 

m, 

N-sec^ 

X, 

m2 

m2 

m 

253.782 

11.4783 

169,328 

12.4663 

321.461 

11.4783 

308.120 

12.4663 

321.461 

11.6053 

308.120 

12.6060 

492.272 

11.6053 

33.012 

12.6060 

492.272 

11.8986 

33.012 

12.6314 

343.931 

11.8986 

190.033 

12.6314 

343.931 

11.9113 

190.033 

13.2280 

410.486 

11.9113 

134.896 

13.2280 

410.486 

11.9177 

134.896 

13.2638 

692.096 

11.9177 

148.099 

13.2638 

692.096 

11.9225 

148.099 

13.3804 

52.869 

11.9225 

163.985 

13.3804 

52.869 

11.9367 

162.378 

13.4899 

127.891 

11.9367 

195.384 

13.4899 

127.891 

11.9621 

193.605 

13.5394 

85.343 

11.9621 

235.532 

13.5394 

85.343 

11.9685 

234.642 

13.6466 

174.651 

11.9685 

126.870 

13.6466 

174.651 

12.0637 

122.230 

13.6857 

164.330 

12.0637 

77.442 

13.6857 

164.330 

12.0955 

76.373 

13.7210 

169.328 

12.0955 

43.540 

13.8480 

169.328 

12.1094 

42.823 

13.9270 

308.120 

12.1094 

58.350 

13.9270 

308.120 

12.2500 

56.744 

14.0208 

95.203 

12.2500 

102.594 

14.0731 

95.203 

12.2618 

102.594 

14.0731 

175.582 

12.2618 

87.432 

14.2227 

175.582 

12.2733 

87.432 

14.2227 

155.339 

12.2733 

78.869 

14.2748 

155.339 

12.3444 

78.869 

14.3560 

169.328 

12.4024 

81.544 

14.3764 

12.4024 

41.038 

14.5659 

12.4151 

43.540 

14.5659 

12.4151 

125.795 

14.6718 

14.6718 

N-sec^ 


133.469 

310.478 

310.478 

356.845 

356.845 
312.443 
312.443 
321.006 
321.006 

32.474 

32.474 

109.846 

109.846 
190.881 
190.881 

37.473 

37.473 

40.148 

38.362 

42.644 

40.679 

35.687 

33.026 

61.025 

60.674 

59.598 

13.383 

10.887 

6.371 

4.826 

4.516 

3.447 

3.213 

11.259 

10.811 

0 


TABLE  2.-  VEHICLE  AXIAL  EXTENSIONAL  COEFFICIENT 
(a)  U.S.  Customary  Units 


in. 

AE, 

Ibf 

0 

295.1  X  106 

1.476 

295.1 

1.476 

280.3 

6.476 

280.3 

6.476 

309.4 

9.076 

309.4 

9.076 

1201.7 

11.376 

1201.7 

11.376 

504.4 

13.976 

504.4 

13.976 

613.6 

15.976 

613.6 

15.976 

699.9 

17.976 

699.9 

17.976 

2350.4 

20.776 

2350.4 

20.776 

841.88 

24.176 

841.88 

24.176 

689.3 

25.376 

689.3 

25.376 

316.7 

25.776 

316.7 

25.776 

576.4 

28.476 

576.4 

28.476 

316.7 

32.376 

316.7 

32.376 

576.4 

37.376 

576.4 

37.376 

316.7 

100.776 

316.7 

100.776 

576.4 

103.776 

576.4 

103.776 

316.7 

161.276 

316.7 

161.276 

576.4 

164.276 

576.4 

164.276 

316.7 

169.376 

316.7 

169.376 

841.9 

171.876 

841.9 

171.876 

316.7 

171.976 

316.7 

171.976 

172.976 

172.976 

176.976 
176.976 
185.356 
185.356 
187.956 
187.956 
193.582 
193.582 
198.5 
198.5 
199.3 
199.3 
203.756 
203.756 
205.916 
205.916 
208.7 
208.7 
209.556 

209.556 
212.256 
212.256 
213.656 
213.656 

215.556 

215.556 

315.556 
315.556 
317.156 

317.156 

320.156 

320.156 

322.156 
322.156 


324.7 

332.2 

335.2 
336.0 
342.0 


1041.3  X  10^ 
1041.3 
2449.2 
2449.2 
39.0 
39.0 
78.0 
78.0 
112.8 
112.8 
135.2 

135.2 

369.2 
369.2 
130.0 
130.0 
118.6 
118.6 
169.0 
169.0 
442.0 
442.0 
280.8 
280.8 
730.1 

730.1 
137.8 
137.8 

70.5 

70.5 

130.0 

130.0 

431.6 

431.6 

213.2 
213.2 

71.42 
98.64 


345.4 

345.4 

346.5 
349.0 
349.0 
350.4 
350.4 

352.3 

452.3 

453.9 
456.0 

456.9 
456.9 
458.8 

458.8 

461.4 
461.4 
461.7 
461.7 

465.6 

467.6 
467.6 

468.9 
468.9 

469.699 
470.011 
470.011 
471.824 
471.824 

489.699 
489.699 
491.035 
491.035 


537.268 

543.268 

549.268 

553.627 

559.627 


571.62 


574.327 

577.629 


199.1  X  10 
612.0 

199.1 

336.7 

730.1 

730.1 
195.0 

70.5 

70.5 

184.5 

529.1 

529.1 

187.8 

187.8 
75.4 
75.4 

707.2 
707.2 

75.4 

75.4 

210.6 

190.9 
452.4 
452.4 
870.7 

870.7 

162.8 

98.5 
33.46 
33.46 

182.16 

182.16 

25.74 


TABLE  2.-  VEHICLE  AXIAL  EXTENSIONAL  COEFFICIENT  -  Concluded 

(b)  SI  Units 


X, 

m 

AE, 

N 

X, 

m 

AE, 

N 

X, 

m 

AE, 

N 

0 

131.266  X  lo"^ 

4.3936 

463.191  X  lo"^ 

8.8011 

88.564  > 

.0375 

131.266 

4.3936 

1089.453 

8.8646 

149.771 

.0375 

124.683 

4.4952 

1089.453 

8.8646 

324.763 

.1645 

124.683 

4.4952 

17.348 

8.9002 

324.763 

.1645 

137.627 

4.7080 

17.348 

8.9002 

86.740 

.2305 

137.627 

4.7080 

34.696 

8.9484 

31.360 

.2305 

534.540 

4.7741 

34.696 

11.4884 

31.360 

.2890 

534.540 

4.7741 

50.176 

11.5291 

82.069 

.2890 

224.367 

4.9170 

50.176 

11.5824 

235.354 

.3550 

224.367 

4.9170 

60.140 

11.6053 

235.354 

.3550 

272.942 

5.0419 

60.140 

11.6053 

83.537 

.4058 

272.942 

5.0419 

164.228 

11.6535 

83.537 

.4058 

311.329 

5.0622 

164.228 

11.6535 

33.539 

.4566 

311.329 

5.0622 

57.827 

11.7196 

33.539 

.4566 

1045.504 

5.1754 

57.827 

11.7196 

314.577 

.5277 

1045.504 

5.1754 

52.756 

11.7272 

314.577 

.5277 

374.485 

5.2303 

52.756 

11.7272 

33.539 

.6141 

374.485 

5.2303 

75.174 

11.8262 

33.539 

.6141 

306.614 

5.3010 

75.174 

11.8770 

93.679 

.6446 

306.614 

5.3010 

196.610 

11.8770 

84.916 

.6446 

140.874 

5.3227 

196.610 

11.9101 

201.236 

.6548 

140.874 

5.3227 

124.905 

11.9101 

201.236 

.6548 

256.394 

5.3913 

124.905 

11.9304 

387.305 

.7233 

256.394 

5.3913 

324.763 

11.9383 

387.305 

.7233 

140.874 

5.4269 

324.763 

11.9383 

72.417 

.8224 

140.874 

5.4269 

61.296 

11.9843 

43.815 

.8224 

256.394 

5.4751 

61.296 

11.9843 

14.884 

.9494 

256.394 

5.4751 

31.360 

12.4384 

14.884 

.9494 

140.874 

8.0151 

31.360 

12.4384 

81.028 

2.5597 

140.874 

8.0151 

57.827 

12.4723 

81.028 

2.5597 

256.394 

8.0558 

57.827 

12.4723 

11.450 

2.6359 

256.394 

8.0558 

191.984 

12.4968 

11.450 

2.6359 

140.874 

8.1320 

191.984 

12.4968 

1.979 

4.0964 

140.874 

8.1320 

94.836 

13.6466 

1.294 

4.0964 

256.394 

8.1828 

94.836 

13.7990 

.850 

4.1726 

256.394 

8.1828 

31.769 

13.9514 

.609 

4.1726 

140.874 

8.2474 

43.877 

14.0621 

.3478 

4.3022 

140.874 

8.4379 

60.865 

14.2145 

.1770 

4.3022 

374.494 

8.5141 

189.315 

14.5191 

.0752 

4.3656 

374.494 

8.5344 

60.865 

14.5879 

2.055 

4.3656 

140.874 

8.6868 

50.087 

14.6718 

0 

4.3682 

4.3682 

140.874 

463.191 

8.7732 

8.7732 

88.564 

272.230 

TABLE  4.-  COMPARISON  OF  SOLUTIONS  TO  EIGENVALUE  PROBLEMS  OBTAINED  FROM 
BOTH  STIFFNESS  AND  FLEXIBILITY  APPROACHES 


7) 

ModX. 

2 

4 

5 

6 

7 

8 

10 

Method^ 

Frequency  ratio 

0.99952 

0.99989 

0.99965 

0.999652 

0.94111 

0.80404 

80.404 

I 

1 

.99952 

.99996 

.99982 

2158.3 

6825.3 

21  580. 

1.5  X  10^ 

II 

.99952 

.99999 

1.0 

1.0 

1.0 

1.0 

1.0 

III 

.99952 

.99999 

1.0 

1.0 

1.0 

1.0 

1.0 

IV 

0.99943 

1.0 

0.99998 

0.99962 

18.001 

I 

0 

.99943 

1.0 

1.0 

1.1202 

.97744 

II 

.99943 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

III 

.99943 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

IV 

0.99979 

0.99999 

0.99999 

0.99992 

1.6042 

13.124 

1310.2 

I 

.99979 

.99999 

.99999 

.99999 

.99988 

.99844 

.99984 

II 

.99979 

.99999 

.99999 

1.0 

1.0 

1.0 

1.0 

m 

.99979 

.99999 

1.0 

1.0 

1.0 

1.0 

1.0 

IV 

1.0 

1.0 

1.0 

1.0 

I 

A 

HHiBM 

1.0 

1.0 

1.0 

II 

*1 

(c) 

(c) 

(c) 

(c) 

(c) 

III 

.99932 

1.0027 

.97831 

(d) 

IV 

Reference  frequency, 

1 

0.65523 

0.65523 

0.65523 

0.65523 

0.65523 

0.65523 

0.65523 

2 

1.3257 

1.3257 

1.3257 

1.3257 

1.3257 

1.3257 

1.3257 

3 

1.8202 

1.8202 

1.8202 

1.8202 

1.8202 

1.8202 

1.8202 

4 

14.142 

141.42 

447.21 

1414.2 

4472.1 

14  142.1 

141  421. 

^  I:  Givens  solution  (ref.  10,  pp.  336-340)  to  stiffness  matrix. 

II:  Iteration  and  sweeping  solution  (ref.  11)  to  stiffness  matrix. 

Ill:  Iteration  and  sweeping  solution  to  flexibility  matrix. 

IV:  Givens  solution  to  flexibility  matrix. 

is  the  calculated  modal  frequencies  of  5  mass  system; 

^ref  reference  modal  frequencies  calculated  from  degenerate  4  mass  and  2  mass  systems. 

^Solution  failed  to  converge. 

^Imaginary  (complex)  frequency. 
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Figure  6.-  Frequency  percent  error  as  a  function  of  total  number  of  lumped  masses. 


Uniform  1:68111  (p  =  O) 


Figure  7.-  Mode-shape  rms  percent  error  as  a  function  of  total  number  of  lumped  masses. 
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”The  aeronautical  and  space  activities  of  the  United  States  shall  be 
conducted  so  as  to  contribute  ...  to  the  expansion  of  human  knowl¬ 
edge  of  phenomena  in  the  atmosphere  and  space.  The  Administration 
shall  provide  for  the  widest  practicable  and  appropriate  dissemination 
of  information  concerning  its  activities  and  the  results  thereof.’' 

— National  Aeronautics  and  Space  Act  of  1958 
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